#include "plot.h"

using std::cout;
using std::cin;
using std::endl;
void output(Para para, int n1)
{
    cout << "test" << endl;
}

//初始化double型数组
void init(double x[], int n)
{
    for (int i = 0; i < n; i++)
    {
        x[i] = 0.0;
    }
}

//初始化complex<double>数组
void init_c(complex<double> x[], int n)
{
    for (int i = 0; i < n; i++)
    {
        x[i] = complex<double>(0.0,0);
    }
}

void output(Para para)
{
    double dx, dz;
    cout << "请输入画图的参数dx(nm)：";
    cin >> dx;
    cout << "请输入画图的参数dz(nm)：";
    cin >> dz;

    double theta[N]; init(theta, N);    //出射角取值数组
    double k[N]; init(k, N);    //kz取值数组
    double epsilon[N]; init(epsilon, N);    //相对介电常数取值数组
    double Rf[N]; init(Rf, N);  //前向反射取值数组
    double Rb[N]; init(Rb, N);  //后向反射取值数组
    double Tf[N]; init(Tf, N);  //前向辐射取值数组
    double Tb[N]; init(Tb, N);  //后向反射取值数组
    complex<double> GRf[N]; init_c(GRf, N); //全局前向反射取值数组
    complex<double> GRb[N]; init_c(GRb, N); //全局后向反射取值数组
    complex<double> fN[N]; init_c(fN, N);   //因子N的取值数组
    complex<double> fO[N]; init_c(fO, N);   //因子O的取值数组
    complex<double> fI[N]; init_c(fI, N);   //因子I的取值数组
    complex<double> sf[N]; init_c(sf, N);   //前向辐射源的取值数组
    complex<double> sb[N]; init_c(sb, N);   //后向辐射源的取值数组
    complex<double> mfN[N]; init_c(mfN, N); //因子N乘积的取值数组
    complex<double> mfO[N]; init_c(mfO, N); //因子O乘积的取值数组
    double **imat = (double **)malloc(sizeof(double *) * NZ);
    for(int i = 0; i < NZ; i++)
    {
        imat[i] = (double *)malloc(sizeof(double) * NX);
    }


    //写入与研究波长相对应的角频率
    double omega = C * 2.0 * M_PI / para.get_wave_l() * 1.0e9;

    //写入粒子速度
    double v = C * sqrt(1.0 / (1.0 + 1.0 / (para.get_p_mon() * 1000.0 * para.get_p_mon() * 1000.0/ (para.get_p_m()*para.get_p_m()))));

    int n = para.get_diel_n() * para.get_r_p_n(); cout << n << endl;

    //写入各域内的相对介电常数数组
    d_perm(para, epsilon, n);

    //数据输出设置
    std::ofstream out;
    out.open("out.dat");
    out << "#域中的前向电场分布" << endl;
    out << "#z" << "\t" << "x" << "\t" << "强度" << endl;

    //开始循环，对于每个出射角theta， 对于每个varphi角
    cout << "进行计算中..............." << endl;
    complex<double> tmp1(0.0), tmp2(0.0); double xx, yy; 
    double z = 0.0;
    double l[N]; para.get_diel_l(l);
    for(int i = 2; i < n+3; i++)
    {
        int index = (i-2) % para.get_diel_n();
        int Nz = l[index]*um/nm/dz + 1;

        //计算theta取值数组
        for(int nk = 0; nk < INT_N1; nk++)
        {
            theta[nk] = asin(1.0*nk/INT_N1);
        }

        for(int nk = 0; nk < INT_N1; nk++)
        {
            kz(epsilon, omega, theta[nk], k, n);
            R1(k, epsilon, Rf, n);
            R2(k, epsilon, Rb, n);
            T1(k, epsilon, Tf, n);
            T2(k, epsilon, Tb, n);
            GR1(para, k, Tf, Rf, Tb, Rb, GRf);
            GR2(para, k, Tf, Rf, Tb, Rb, GRb);
            facN(para, k, Tf, Rf, GRb, fN);
            facO(para, k, Tb, Rb, GRf, fO);
            facI(para, n, GRf, GRb, k, fI);
            source1(omega, v, para, k, epsilon, sf);
            source2(omega, v, para, k, epsilon, sb);
            tmp1 = omega/C*sin(theta[nk]) * Ezf(i, n, sf, sb, GRf, GRb, fN, fO, fI); //cout << i << endl;

            
            for(int nz = 0; nz < Nz; nz++)
            {
                for(int nx = 0; nx < NX; nx++)
                {
                    for(int nphi = 0; nphi < INT_N2; nphi++)
                    {
                        xx = nz*dz*nm*k[i-1]+omega/C*sin(theta[nk])*cos(1.0*nphi/INT_N2*2*M_PI)*nx*dx*nm;
                        tmp2=complex<double>(cos(xx), sin(xx));
                        imat[nz][nx] = (tmp1 * tmp2).real() + imat[nz][nx];
                    }
                }
            }   
        }
        
        for(int nz = 0; nz < Nz; nz++)
        {
            for(int nx = 0; nx < NX; nx++)
            {
                out << nz*dz+z << "\t" << nx*dx << "\t" << imat[nz][nx]*2*M_PI/INT_N2*(omega/C)/INT_N1*1.0e-168 << endl;
                imat[nz][nx] = 0.0;
            }
            out << endl;
        }
        
        
        z += l[index];
    }

    free(imat);

    out.close();
    cout << endl << "数据写入完成********" << endl;
}
